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REMARKS /ARGUMENTS 

Claims 1-10 are in the application for 
consideration . 

1. The examiner has made an objection to Fig. 7 of 
the drawings and has stated that a legend such as "PRIOR 
ART" should be used to designate it '"...because only 
that which is old is illustrated." 

Fig. 7 of the present application, which 
illustrates an example of a configuration of discrete 
image sensing elements which may be employed in 
accordance with certain elements of the present 
invention, has been incorporated in the application to 
assist those skilled in the art to better understand the 
claimed invention , 

Enclosed is a proposed corrected copy of Fig. 7 
showing a proposed change thereto. Applicants propose to 
insert the term ''PRIOR ART", as shown in red, on the 
sheet as suggested by the examiner. This proposed change 
overcomes the objection made to Fig. 7. 

2. Claims 1-10 have been rejected under 35 U.S.C. 
§ 103(a) as being unpatentable over U.S. Patent 
6,181,376 Bl ("Rashkovskiy") . 

Applicants traverse this ground of rejection. The 
reference does not teach or in any way suggest within 
the meaning of Section 103 a method or apparatus for 
electronically capturing and processing image 
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information which includes the method steps of 
applicants or, further, which involves the application 
of two one-dimensional non-linear interpolation 
processes . 

The method and apparatus of applicants involve 
the application of two one-dimensional non-linear 
interpolation processes. Generally speaking, applicants' 
method essentially "decouples" the conventional two- 
dimensional process into two one-dimensional processes. 
The first one-dimensional color recovery process 
generates intermediate second color image data from the 
sampled first color image data and the second one- 
dimensional color recovery process generates the desired 
third color image data from the second color image data. 

As described in detail in the application 
(see, for example, page 8, line 28 to page 10, line 14) 
a two-dimensional array of discrete image sensing 
elements, each of which is specifically responsive to 
one of at least three predetermined colors (e.g., red, 
green and blue) is exposed to image information-bearing 
illumination to obtain a collection of each electronic 
information signal received from each discrete element. 
The collection of signals forms the raw unprocessed one- 
color image data from which fully recovered third color 
image data can be derived. 

The first sequence of steps of applicants 
towards deriving the fully-recovered third color image 
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data is directed to first recovering missing color 
information along a first dimension (e.g., along rows of 
the array) by interpolating the first color image data 
along the first dimension to provide first-interpolated 
color data for each of the discrete elements [see (c) (i) 
of claim 1]. This interpolation step may be carried out 
with any type of filter. 

A difference channel is then formed between the 
first color image data and the first interpolated color 
data [see (c) (ii) of claim 1]. A one-dimensional non- 
linear filter is then applied on the difference channel 
whereby the first-recovered image data are obtained as a 
combination of the first color image data and the 
filtered first difference channel [see (c) (iii)of claim 
1] . The first color data are then combined with the 
first-recovered color data to obtain the second color 
data [see (c) (iv) of claim 1]. 

This sequence of method steps is basically 
repeated using the second color data to derive fully- 
recovered third color image data from the two color 
image data by recovering missing color information along 
a second dimension, e.g., along columns of the array. 
The second recovered image data are obtained by 
interpolating the second color image data along the 
second dimension to provide second interpolated data for 
each of the discrete elements, forming a difference 
channel between the second color image data and the 
second interpolated data and then applying a one 
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dimensional non-linear filter on the difference channel 
and combining with the second color data. The third 
color image data comprise a combination of the second 
recovered image data and the second image data. 

The method of Rashkovskiy does not teach or suggest 
the method steps of applicants and nowhere teaches or 
suggests using a non-linear filter in a method to 
recover missing image data. 

The reference does teach a method to generate 
missing color values for pixels in a color filter array 
created by a digital camera but does not teach or 
suggest the method and apparatus of applicants. More 
specifically, the reference does not perform 
steps(c) (ii), (c ) (iii) and (c) (iv) or (d) (ii), (d) (iii) 
or (d) (iv) of applicants' method or use a non-linear 
filter in any step. 

Generally speaking, the first step of the 
Rashkovskiy method is to find the green values at all 
red and blue sites where green was not measured (see, 
for example. Fig. 4a and Fig. 4b). Interpolation is 
performed to find the green values at the diagonal 
corners using a cubic B-spline filter followed by 
another cubic B-spline interpolation filter in the 
orthogonal diagonal direction (column 4, lines 23 - 46). 

The examiner' s statement that a cubic B-spline 
filter is considered to be a non-linear filter is not 
correct; a cubic B-spline filter is a linear filter. 
This is shown unequivocally by Splines, A Perfect Fit 
for Signal and Image Processing, M. Unser, IEEE Signal 
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Processing Magazine, November 1999, pages 22-38. At 

page 23, following equation (3), the article states 

The B-splines of degrees 0 to 3 are shown in 
Fig. 1. Since the B-spline model (1) is linear 

Although a cubic B-spline is not specifically mentioned 

by name in the above quoted statement, a B-spline of 

degree 3 is in fact a cubic B-spline. This is seen 

clearly in equation (6) on page 24 of the article, which 

is described as providing a ^'closed-form representation 

of the cubic B-spline . . .which is often used for 

performing high quality interpolation." The superscript 

3 used for the B-spline in equation (6) indicates that a 

degree 3 B-spline is a cubic B-spline. A copy of this 

article is enclosed for the convenience of the examiner. 

It should also be noted here that the Rashkovskiy 

reference cites an article by a group of authors 

including Unser for cubic B-spline filters (see column 
« 

4, lines 23 -28) . 

Further at column 5, lines 59 - 65 and column 6, 
lines 4-7 of the reference, specific and general 
formulae are provided for the B-spline interpolation 
filters. The formulae are for linear filters as is 
apparent to those skilled in the art since the 
interpolated values are obtained as a weighted linear 
sum of the neighboring pixel values. 

Having obtained the missing green data the 
method of Rashkovskiy does not form second color data as 
is obtained according to the method of applicants as 
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described in detail above, but rather goes on directly 
to obtain missing red values from data at the corners, 
again by linear interpolation. Broadly stated, the 
reference teaches going directly from the first linear 
interpolation step to a second linear interpolation 
step. The method does not form any intermediate color 
data as is the case in the method of applicants. 

Furthermore, all formulae specified in the 
Rashkovskiy reference for recovering the missing red and 
blue colors, subsequent to the cubic B-spline 
interpolation of green, involve only linear combinations 
of neighboring pixel values. These operations are 
therefore linear. Since cubic B-spline interpolation is 
also linear, the Rashkovskiy reference does not involve 
any non-linear filtering operation to recover the 
missing colors. 

It is apparent from the foregoing remarks that the 
reference does not teach or suggest applicants' claimed 
method and apparatus. For all the foregoing reasons, 
claims 1-10 are clearly patentably distinguishable 
over the disclosure of Rashkovskiy. Those skilled in the 
art would find no teaching or suggestion in this 
reference which would place them in possession of 
applicants' claimed method and apparatus. 

Reconsideration of this ground of rejection and 
withdrawal thereof are respectfully requested. 

In summary, it has been shown that the claims are 
directed to wholly novel and unobvious subject matter. 
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Reconsideration of the application and allowance of the 
claims are respectfully requested. 
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A Perfect Fit for Signal and Image Processing 



[ hiding a general mechanism for switching be- 
tween the continuous and discrete signal do- 
mains is one of the fundamental issues in signal 
processing. It is a question that arises naturally 
during the acquisition process where an analog signal or 
image is to be converted into a se- 
quence of numbers (discrete repre- 
sentation) . Coiiverscly, the need for a 
conrinuous signal representation co- 
mes up every time one wishes to im- 
plement numerically an operator that is initially defined 
in the continuous domain. Typical examples in image 
processing are the detection of edges through the com- 
putation of gradients (spatial derivatives), and geomet- 
ric transformations such as rotations and scaling 
(interpolation) . 



Michael Onsei 



The textbook approach to those problems is provided 
by Shannon's sampling theor)s which describes an equiv- 
alence bet^^'een a band-limited function and its equidis- 
tant samples taken at a frequency that is superior or equal 
to the Nyquist rate [76] . Even though this theory has had 

an enormous impact on the field, it 

still has a number of problems. 
First, it relies on the use of ideal fil- 
ters, which are devices not com- 
monly found in nature. Second, the 
band-limited hypothesis is in contradiction widi the idea 
of a finite (or finite duration) signal. Third, the 
bandlimiting operation tends to generate Gibbs oscilla- 
tions, which can be visually disturbing in images, Finally, 
the underlying cardinal basis function [sine (a?)] lias a very 
slow decay, which makes computations in the signal do- 
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main very inefficient. While the £ii*st two problems can be 
dealt with by using approximations and introdxicing con- 
cepts such as an esscnUat bandwidth and an essential tim 
duration [78], there is no way to address die last two is- 
sues odier than changing basis functions. 

Our purpose here will be to provide arguments in fa- 
vor of an alternative approach that uses splines, which is 
equally justifiable on a theoretical basis, and which offers 
many practical advantages, To reassure the reader who 
may be afraid to enter new territory, we must emphasize 
that we are not losing anything because we will retain die 
traditional dieory as a partiailar case (i.e., a spline of infi- 
nite degree). The basic computational tools will also be 
fiamiliar to a signal processing audience (filters and recur- 
sive algoritlims)> even though their use in the present 
context is less conventional. In the course of the presenta- 
tion, we will also bring out the connection with the 
mulriresolution theory of the wavelet transform. 

Interestingly, splines are slightly older than Shannon's 
sampling theorJ^ They were first described in 1946 [70]. 
In his landmark paper, Schocnberg laid the mathematical 
foundations for the subject; he showed how one could 
use splines to interpolate equally spaced samples of a 
function, He also introduced the B-splines, the basic at- 
oms by which polyiiomial splines ai'c constructed. De- 
spite this early start, the subject of splines then lay more or 
less dormant during tiie 1950s, while signal processing 
developed at a rapid pace within Shannon's elegant 
framework of band-limited functions. Splines only really 
took off in the early 1960s when mathematicians realized 
that diesc functions could model the physical process of 
drawing a smoo± curve (minimum curvature property). 
This a-eated an intense interest in the subject and the ap- 
plications soon followed in approximation theory [24], 
[74], numerical analysis [64], and various other branches 
of applied mathematics [3]. With the advent of digital 
computers, splines caught the interest of engineers and 
had a tremendous impact on computer-aided design 
[29], [45] and computer graphics [10]. Howe\^er, there 
was littic crossover to signal processing, perhaps because 
researchers in this field had become so accustomed to 
thinking in terms of band-limited functions. Recentiy, 
thanks in part to a nevv (non-band-limited) way of think- 
ing brought forth by wavelet theory [51], die situation 
has changed significantiy. 

This article attempts to fullfill three goals, The first is 
to provide a tutorial on splines that is geared to a signal 
processing audience. The second is tx> gather all their im- 
portant properties and provide an overview of the mathe- 
matical and computational tools available; i.e., a road 
map for the practitioner with references to the appropri- 
ate literature. The third goal is to give a review of the pri- 
mary applications of splines in signal and image 
processing; most of those are discussed in the final part of 
the article. 



Spline Interpolation 
Polynomial Splines 

Splines are pieccwise polynomials with pieces that are 
smoothly connected together. The joining points of the 
polynomials are c^kdhtots. For a spline of degree n, eacli 
segment is a polynomial of degree », which would sug- 
gest that we need « -I- 1 coefficients to describe each piece. 
However, diere is an additional smoothness constraint 
tiiat imposes the contmuity of the spline and its deriva- 
tives up to order (w-1) at the knots, so that, efTcctively, 
there is only one degree of freedom per segment. Here, 
we will only consider splines with uniform knots and unit 
spacing. The rcmarlcable result, due to Schoenberg [70], 
is that these splines are uniquely characterized in temis of 
a B-spline expansion 



(1) 



which involves the integer shifb of the central B-spline of 
degree w denoted by (J "(.«); the parameter^ of the model 
are the B-spline coefficients c(k). B-splincs, defined below, 
are sjTnmetrical, beJl-shaped functions consoncted from 
the (»-|-l)-fold convolution of a rectangular pulse p**: 



1, 
1 

2 > 

0, 



otherwise 



p'*(x)=P''*P'>»-*p°(*). 



(2) 



(3) 



{ n+lj times 



The B-splincs of degrees 0 to 3 are shown in Fig. 1. Since 
the B-spline model ( 1 ) is linear, studying the properties of 
the basic atoms can tell us a lot about splines in general 
(cf, Box 1). Thanks to this representation, each spline is 
unambiguously cliaracterizcd by its sequence of B-spline 
coefficients c{k)^ which has the convenient structure of a 
discrete signal, even though the underlying model is con- 
tinuous (discrete/continuous representation). 




A l,The centered B-spfines of degree 0 to 3. 
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B -splines are very easy to manipulate. For instance, we 
can obtain derivatives through the following formula 



fbc 



(4) 



which reduces the degree by one. Similarly, we compute 
the integral as 



(5) 



Once we know the effect of linear operators sucli as (4) or 
(5) on the basis functions, it is a trivial matter to apply 
them to any spline via tlie B-spline representation (1). 

Within the family of polynomial splines, cubic splines 
tend to be the most popular in applications — perhaps due 
to their minimum curvatm-e propert)', which is discussed 
in Variational Properties." Using (2), we obtain the fol- 
lomng closed-form representation of the cubic B-spline 



l<^<2 
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0, 



(6) 

which is often used for performing high-quality inter- 
polation. 



B-'SpNne Interpolation via Digital Filtering 

From what has been said so far, it appears that most of the 
work consists in determining the B-spline model of a given 
input signal We now consider the spline interpolation 
problan where the coefHcients are determined such that 
the function goes dirough the data points exacdy (cf Fig. 
2). For splines of degree 0 (piecewise constant) aixi splines 
of degree 1 (piecewise linear), this is a trivia! matter be- 
cause the B-spline coefficients are identical to the signal 
samples: c(k)—s{k). For higher-degree splines, however, 
the situation is more complex. Traditionally, the B -spline 
interpolation problem has been approached using a matrk 
framework — setting up a band -diagonal system of equa- 
tions, whicli is then solved using standard numerical tech- 
niqucs (forward/backward substitution or LU 
decomposition) [25], [64]. In the early 1990s, it was rec- 
ognized that this problem (as well as many other related 
ones) could also be approaclied using simpler digital filter- 
ing tecliniques [33], [93], [96], [97]. 

To derive this tj^pe of signal processing algorithm, we 
need to introduce the discrete B-spline kernel bl , which 
is obtained by sampling tlie B -spline of degree n ex- 
panded by a factor of m\ 



(11) 



;Bpxi.B- 





iSSittSiSiiiiiliSiM 

^^bfabtdrs^(tune^5hifu)^and dp-.;: 




24 



IEEE SIGNAL PROCESSING MAGAZINE 



NOVEMBER 1999 



















.3'- 
















. 2. 
















■ V 




































' Cubic Spline . 
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^ Cubic BrSpnne Baists Functions 





A 2, Example of a cubic spline signal that is represented as a lin- 
ear combination of shifted cubic B-splines. Given the signal's 
samples, the basic problem is to determine the appropriate co- 
efficients in (i). 

N0W5 given tlie signal samples s{k)^ we want to deter- 
mine the coefficients c{k) of the B-spline model (1) such 
that we have a perfect fit at the integers; i.e. » Vife eZ, 

Using the discrete B-spiines, this constraint can be rewrit- 
ten in the form of a convolution 

s[k)^{b!:^c){k\ (12) 

Defining the inverse convolution operator 

the solution is found by inverse filtering (cf [97]) 

c{k)^{bn-'^s{ky (13) 

Since hi is symmetric PIR (finite impulse response), the 
so-called direct B-spline filter {b\ is an all-pole system 
tliat can be implemented very efficiently using a cascade 
of first-order causal and anti-causal recursive filters [93], 
[96] , Tliis algorithm is stable numerically and is faster and 
easier to implement than any other numerical tecliniquc. 



The explicit procedure for the cubic spline case is de- 
scribed in Box 2. 



Cardinal Splines 

To bring out the connection between the spline interpo- 
lation process and the traditional approach for 
band-limited functions, it is helpful to introduce the car- 
dinal spline basis functions that are the spline analogs of 
the sine function. Combining (1) and (13), we have 

xw=g(w)-'*0(*)P"(*-^) 
■=S^(*)SW)"'(')P'(*-'-*) 

*e2 (14) 

where we have identified the cardinal spline of degree n : 



(IS) 



Thus, (14) provides a spline interpolation formula diat uses 
die signal values as coefficients. The formula works because 
Tj " (x) has the same interpolauon property as the sine func- 
tion; it vanishes for all integers except at the or^in, where it 
takes the value one, The cardinal spline represents the im- 
pulse response of the corresponding spline interpolator, 
is^ote that, for n ^ 2, these fiincdons are no longer compacdy 
supported; however, they decay exponentially fast. We can 
also express (15) in the Fourier domain, which yields the 
frequency response of the spline interpolator of degree n 



H 



1 



(16) 



The cardinal cubic spline is shown in Fig, 4 and appears to 
be quite similar to the sine function. In fact, it has been 
shown that ri " (x) converges to sinc(x) as w goes to infinity 
[7]. It is a rather strong type of conveigeiice (I^-norm) 
that holds in both time and frequency domains (cf. Fig. 5) . 

Note that the correspondence between splines of infi- 
nite order and band-limited functions was Icnown to 
Schoenberg and his successors [27], [73]. However, 
these mathematical results did not reach tlic signal pro- 
cessing community until recently [7], mainly due to sub- 
stantial differences in context and terminology. 
Approximation theorists t)picaily speak of "entire func- 
tions of exponential type" when they refer to 
band -limited functions, The recent cross-fertilization that 
has occurred has been quite fruitful and there have been 
benefits on both sides. For instance, the idea of using an 
anti-aliasing filter in Shannon's sampling tlieory has sug- 
gested similar solutions for splines; these are discussed in 
the next section. 

We should emphasize that the primary usefulness of 
the cardinal splines is conceptual. They provide us with a 
better understanding of the algorithm. From a practical 
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point of view, however, it is much more efficient to work 
with the B -spline representation, at least when we are per- 
fonning interpolation. The reason for this is that, in most 
applications J it is die re -sampling part (evaluation of the 
expansion formula (1) or (17)) that is by far tiic most 
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A 5. Recursive causai and anthcausal filters for cubic spline inter- 
polation. 
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A 4. Cardinal (or fundamental) cubic spline. 




A 5. Frequency response of the spline interpolators of degree 
n^l and 3,Asn increases, the spline interpolators fend to the 
ideal low-pass filter (dotted line), 

costly step. Accordingly, we have the advantage of using 
the shortest possible basis functions (i,e., B-splines) such 
that the number of terms that contribute for a given ^ is 
minimized. This is precisely why splines are so much 
more computationally efficient than the traditional 
sinc-based approach, Because sinc(A:) decays lilce l/\x\y 
computing a signal value at a particular non-integer loca- 
tion with an error of less than 1% will require of the order 
of 100 operations in each direction, while B-splines pro- 
vide an exact computation with just a few non-zero terms 
(«+ 1 to be precise) . (In ^-dimensions, the complexit)'^ of 
an interpolation algorithm that uses separable basis func- 
tions increases with die power of^. For diis reason, virtu- 
ally no one uses sine interpolation for images, not to 
mention volumes.) An illustration of how these ideas can 
apply for the geometric transformation of images is given 
in Box 3. When compared to any other type of 
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interpolator, B-spUnes offer the best performance for the 
least complexity (we discuss this issue in more depth later 
in the article). 

It is worth mentioning that niany authors in image 
processing leave out the essential prcfiltering step in Box 
3» This has a catastrophic effect on performance and per- 
petuates the incorrea belief that high-order B-spline in- 
terpolation results in increased image blurring. By 
looking at the frequency responses in Fig. 5, the reader 
will be convinced ^at this cannot be the case. In fact, us- 
ing a higli-degrec spline interpolator (typ*, «=5) is a 
good, stable coniputadonal way of approximating the 
ideal sine interpolator. For a given computational bud- 
get, it is usvially superior to using a windowed sine func- 
tion, especially in high dimensions [84]. The spline 
approach wins because it has a high order of approxima- 
tion (cf, "Controlling the Approximation EiTor") and an 
eflfective, sinc-like impulse response that, is infinite 
(tlianks to the IIR prefiltering step). 



Spline Sampling Theory 

Most of die devebpments in this area are relatively recent 
and have greatly benefited from the analogy of the tradi- 
tional approach dictated by Shannon's sampling theory 
which recommends the use of an anti-aliasing filter when 
the input signal is not band-limited [38], [95] . The concepts 
are best explained from die general perspective of Hilbert 
spaces [6] , For convenience, we will use a slightly more gen- 
eral spline generating duiction, whicli we represent as 

9W = I^WP"(A^-^), (18) 

NOVEMBER 1999 



with the important restriction that the sequence p is 
such that the integer translations of q) form a basis of 
our basic spline space, (The necessary and sufficient 
condition for having a Ricsz basis is that 
0<|?(^^"')|<-H», where P{c^) is the Fourier transform 
of ±e sequence^(j&) [61). The two. special cases that we 
have in mind are the B-splines, with p(:k)=^[k] (the 
Kroncckcr delta), iand the cardinal splines , with 
p = (bi)~^. Since we are also interested in varying the 
sampling step, we defme the spline space of degree w 
with step size r by rescaling the. basic model in (1) 

I ". J ■• (19) 

These splines with step size T are formed by taking linear 
combinations of the generalized spline basis functions 
rescaicd by a factor T and spaced accordingly. As before, 
±ere is exactiy one coefficient c{k) per knot or sampling 
point. The condition c(j&) elj means that we are restricting 
ourselves to linear combinations with a finite energy. In 
this way, we are ensuring diat Sr Is a : well-defined 
subspace of Lj, the space of all finite energy funaions. 
Note that the space is considerably larger than Bz = , 
die traditional subspace of band-limited functions con- 
sidered in signal processing. To use an analogy, L3 is to Bn 
(or S{) as R (die real numbers) is to Z (the integers). 

Spline Sampling via an Appropriate Pref liter 

Now, we are interested in approximating an arbitrary sig- 
nal s{x) by a spline 5 ^ eS? . As a measxire of error j we use 
the Xj-norm \\s - || , which is induced by the X^^i^iner 
product: 
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s{x) 





Sampling 


1 ( lUI^MPI^ 









► 



A & Least-squares spi'me approximation. The analog input signal 
s(x) is prefiftered with ©(-x) and sampled thereafter to yield 
the generalized spline coefficients c(k)^(^ s(x),^(x-k) ). The 
sampling is modeled by a multiplication with a sequence of 
Dirac deltas. The spline approximation s(x) = J^^^^^ c(/r)tp(x-fr) 
Is then obtained by post-filtering with <p(x). 
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A 7. Frequency response of the optimal spline prefilters of degree 
n=l and J. As n increases^ the optimal spline prefilters tend to 
the ideal low-pass filter (dotted line), 



(20) 



According to tiiis criterion, the minimum error approxi- 
mation of s(x) 6 13 in 5; Is given by its orthogonal pro- 
jection onto S^. Using the property that the 
corresponding approximation error s(x) - Sr (x) must be 
orthorgonal toS ", it is not too difficult to show that die 
coefficient of the best approximation (least-squares solu- 
tion) are given by (cf. [95]) 



(21) 



where cp(A') eS{ is die dual of (p(x), in the sease that 



(bi-orthogonality condition). While this may sound 
rather abstract, diere is a simple prefiltering and sampling 
interpretation of (2 1 ). The corresponding block diagram 



for the simplified case of a unit sampling step T= 1 is 
shown in Fig* 6. 

This is similar to the conventional sampling procedure 
dictated by Shannon's clieory except that the optimal 
prefilter, which is the time-reversed version of <p(jc), is not 
necessarily ideal. In the particular case of the cardinal rep- 
rcscntadon where the spline coefficients are the signal 
samples (i.e., p = {bi ), the transfer function of die opti- 
mal prcfilter is given by (cf. [95]) 



/f(a)) = 



^ sin(Q)/2) y^ B^je^'') 



{22} 



where B{(e^'^) is the Fourier transform of a discrete 
B-splinc of degree j& [cf. (11)]. 

The frequency responses of the optimal prefiltcr for 
the cardinal spline representations of increasing dcgi-ees 
are shown in Fig. 7. The low-pass character of die re- 
sponse suggests that the prefilter ^{x) has a role analo- 
gous to that of the anti-aliasing filter required in 
conventional sampling theory. In fact, as tj;ic order of the 
spline goes to infmit)^ both H " ico) and H " (co) [cf (16) 
and (22)] converge to the ideal low-pass filter (dotted 
lines in Fig. 7) [95], which is consistent with the fact that 
a band-limited signal can also be viewed as a spline of infi- 
nite degree (i.e., Bn =S^ ). 



Controlling the Approximation Error 

We have just seen that there is no fundamental difference 
between the process of performing a least-squares spline 
approximation of a signal and obtaining its band-limited 
representation using the standard sampling procedure 
dictated by Shannon's theory. The only difference is in 
the choice of die appropriate analog prefilter. So far so 
good, but how should we choose the sampling step T> Is 
there any equivalent of the sampling theorem that tells us 
tiiat the signal can be reconstructed exacdy if it is sampled 
at a frequencj^ 1/T that is at least twice the Nyquist rate 
cOtnax /(27i)? In principle, one should expect a similar re- 
sult, at least for higher-order splines. 

Because we are performing an orthogonal projeaion, 
the approximation error wiQ be generally non-zero unless 
the signal is already included in our approximation space. 
However, we can hope to control tliis error by choosing a 
sampling step T that Is sufficiendy small. To'analyzc this 
situation, which is more complicated than in the tradi- 
tional band-limited case, we turn to approximation the- 
ory. A fundamental result is that the rate of decay L of the 
error as a flmction of Tdepends on the ability' of the rep- 
resentation to reproduce polynomials of degree n=i-l. 
The approximation error also depends on die bandwiddi 
of the signal. The relevant measure in this context is 



" [2ni 



f(co)|^ day 



(23) 
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where s{(£i) denotes the Fourier transform of this is 
nothing but tlie norm of the Ldi derivative of s. The key 
result from the Strang-Fbt dieory of approximation is the 
following error bound (cf, [40, SO]): 



£'\'G)) = l-lf''(CO)H''((fl). 



(24) 



where Pr j is the least-squai'cs spline approximation ofs at 
sampling step Tand Q is a known constant. TF"/ denotes 
the space of functions that are I times differentiable in the 
X^, or finite-energy sense. In other words, the error will 
decay like 0(r^\ where the order I=»+l is one more 
than the degree ». Spline interpolation gives the same rate 
of decay as the least-squares approximation (21), but 
wida a larger leading constant [103]. 

Recently, it has become possible to determine the ap- 
proximation error much more precisely by simply inte- 
grating the whole spectrum of the function to 
approximate against a frequency kernel £ * (co) [15]. The 
justification for this procedure is the error fonnula 



ll^-Pr^ll = 



2n 



(25) 



where |y| is bounded by some known constant [13]. The 
second term in (25) is a correction tiiat may take positive 
or negative values. It is zero for band-limited fiinctions 
and very small otherwise, provided thAts{x) is sufficiently 
smooth {s{x)eW2mth r large). Moreover, die second 
tei'm cancels out if we take the average approximation er- 
ror over all possible shifts of the input fimction; this is a 
reasonable thing to do since tiie sampling phase is usually 
arbitrarj\ Thus, the first term on the right-hand side in 
(25) provides a very accurate prediction of the error, 
which can be the basis for a quantitative Fourier domain 
evaluation [15]. The error kernel for a least-squares spline 
approximation of degree n is 




k. a Frequency plot of the error kernels for the feast-squares 
spline approximations of degree n«0, /,2, J. Below the Nyquist 
frequency ()>=^ the kernels (and therefore the errors) tend to 
be smaller for splines of higher degree. 



(26) 



where H " (co) and H " (co) are the spline filters defined by 
(16) and (22), respectively. The main point is that tiic 
study of these kernels gives us a very direct way to assess 
the performance of the various types of spline approxima- 
tions. Tills approach is simple, inuiitive, and yet powerful 
enough to recover all classical results and X^-bounds in 
approximation theory [e.g. (24)], We have plotted tiie er- 
ror kernels for «=0, to 3 in Fig, 8. This graph clearly 
shows that for signals that are predominantiy low-pass 
(i.e., witia a frequency content within the Nyquist band) , 
the error tends to be smaller for higher-order splines. 

llie order property (24) is a direct consequence of the 
degree of flamess of the kernel around the origin. Spe- 
cifically, for anith order spline, E " (co) has 2L - 1 vanish- 
ing derivatives atO)=0. This implies that 

£"(rco) = (C2.T^-co'')' +O(r^^*^co^^*^) 

as (0 ^ 0, which explains die 0(T ^ ) behavior of the error 
described by (24) (for more details, refer to [15]). 

As the graph in Fig. 8 also suggests, the error ap- 
proaclies that of a band-limited approximation as the or- 
der of the spline increases, which again reinforces the 
analogy mxh Shamion's sampling theorem. In the limit as 
w-^+oo, the product ii"(co)H'(cD) tencis to the ideal 
low-pass filter, so that we end up with an error entirely 
due to the out-of-band frequency content of die signal. 
The implication is that higher-order splines will usually 
produce better approximations in the -norm, although 
this may occur at the expense of ringing artifacts as the 
model gets closer to being band-limited. 

Multiresolution Spline Processing 

Consider a spline with knots at the integers and dilate it 
by an integer factor m. The resulting enlarged function is 
clearly piecewise polynomial in each unit interval, which 
means that it is also a spline with respect to the initial inte- 
ger grid. This simple observation is the key to the 
multiresolution properties of splines, which maizes them 
perfect candidates for the construction of wavelets and 
P)Tamids. Flere. we will emphasize die special two-scale 
relation for splines, and the construction of pyramids (we 
touch on the subject of wavelets only briefly here). We 
will also briefly make the connection bcDA^een splines and 
wavelets. 



ffl*5ca/e Relation 

For the above scale -invariance argument to hold, we need 
the spline knots to be positioned on the integers. To sim- 
plify the discussion, we will momentarily consider the 
shifted causal B -splines 



cp"(^) = P^(^v-^), 



(27) 
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which have the required property. Similar to the centered 
B-splincs (3), these can also be constructed from die 
(« + l)-fold convolution of cp^. the indicator function in 
the unit interval. Qearly, (^^(x/m\ which is one for 
A' e[0,w), and zero otherwise, can be written as 

bTo (28) 

where ^w(^)jis the filter whose 2;-transform is 

■^"*(^) = Xm ('^^^^'^^^^ pulse of size m). By 
convolving tfiis equation with itself (w + l)-times and per- 
forming the appropriate normalizationj one fmds that 

r(^hn)^J,hZ{ky^''{x-k\ (29) 

where 

This is a uvo -scale equation, which indicates that a 
B-spline of degree n dilated by w can be expressed as a lin- 
ear combination of B-splines. Widi the appropriate phase 
shift, this result also carries over for centered B-splines of 
degree n odd; an alternarive proof is given in [101] . There 
are two remarkable facts connected to the above result. 
First, the two-scale equation (29) holds for any integers 
— not just powers of two, as encountered in the 
muldresolution theory' of the wavelet transform [48], 
[81]5 [109]. vSecond, the refineihcnt filter is simply the 
(«+l)-foid convolution of the discrete rectangular im- 
pulse of width m\ this can be the basis for some very fast 
algorithms [101]. In the standard case where w=2, 
H2 (s) is die celebrated binomial filter that plays a crucial 
role in die theory of the wavelet transform [81] . The filter 
coefficients appear in the Pascal triangle represented on 
the first page of this article, The two-S(ic relation is illus- 
trated in Fig. 9 for the case of the centered B-spline of de- 
gree 1; this corresponds to the third line of PascaFs 
ti'iangle. 

Spline Pyramids 

For constructing multiscale representations of signals, or 
p}Tamids, one usually considers scaling factors that arc 
powers of tw^o. The implication of the rwo-scale relation 
for ?»=2 is tiiat the spline subspaces 5i',with w=2',are 
nested: dSj" d-oSj",--. 





,1 

I ■ ■ 




■ , ■ 1 


-2 


m=2 


•*-2 
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A 9. Illustration of the two-scale relation for the linear B-spline. 



Let P^iS = Si denote the minimum error approxima- 
tion of some continuously defmed signal s(x)eL2 at the 
scale m =2' . We choose to represent it by the following 
expansion 

P,^s^^c,.m{^/2*-k\ (31) 

[or the (p(^/2' -k)& are the spline basis functions at the 
scale (B-spline or odiers); diey are enlarged by a 

factor of 2 * and spaced accordingly] . The expansion coef- 
ficients p^i (k) are defined, at least formally, through the 
inner product (21), The interesting implication of the 
splme nestedness property is that the coefficients c^i (k) 
can be computed itenrivcly in a very simple fashion using 
a combination of discrete prefiltering and 
down-sampling operations. The key observation is that 
we can obtain P^^ ^ = if we simply reapproximate at 
die next finer scale (i.e., P^iS-F^iSi.^). Thus, we may 
compute the expansion coefficient as [cf. (21)] 

^2' (^) = Tr(X^2-<P(^/2"' -l)M^/2' -k)\ . 

2 I (32) 

Using die two-scale relation to precompute the sequence 
of inner products, 

(33) 

it is not difficult to show that the c^ -, {k) are evaluated by 
simple prefiltering widi/; and down-sampling by a factor 
of two: 

c,^{k)^hc,.,){2k\ 

There is also a complementary "interpolation" filter diat 
allows the extrapolation of a coarser resolution to the next 
finer one. An example of such a pyramid is shown in Fig. 
10, where we have used a cardinal representation; in 
other words, we are displaying the samples of the under- 
lying spline images. The corresponding 2D spline model 
is separable, and the procedure is implemented by succes- 
sive 1 D filtering and decimation of the raws and columns 
of die image. The error arrays on the right are obtained by 
subtracting the next coarser approximation from the cur- 
rent spline approximation; it displays the loss of informa- 
tion introduced by image reduction. Specific filter 
formulas can be found in [99]; die filter coefficients and. 
ID approximation routines in the C language can also be 
obtained from the audior on request. 

Instead of minimizing the continuous Lj-error, it is also 
possible to construa splbie pyramids diat are optimal in 
the discrete l^-notm [8]; practically, this amounts to a 
small modification of the reduction filter A [96]. This kind 
of algorithm provides an efficient fdter-based implementa- 
tion of the technique known as spU>ie regression in statistics 
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[2, 31, 113]. Most of the spliiie pyramids use symmetric 
filters that arc centered on the origin (in fact, these arc 
based on the centered B-splines ratlicr that die causal ones 
that have been used here to simplify the argument). Re- 
cently, there has been some incentive for designing centered 
pyramids where die coarsei' nodes arc positioned in die 
center of the finer ones [17]. Hicse last structures are espe- 
cially useful when dealing with image labels. It has also 
be^n shown that shifting the spline luiors between levels 
can improve energy compaction [55]. 

Spline Wavelets 

Tlie L2-splinc pyramid that has been described above has 
all the required properties for a multircsolurion analysis 
oflj in the sense defined by Mallat [50], [51]. In particai- 
lar, the error bound (24) guarantees that we can approxi- 
mate any L^-function as closely as \vc wish by letting the 
scale go to zero. In the wavelet terminology, the 
multircsolurion analysis is dense hiLj [51], HcncCj there 
is no maj{>r difficiilty in constructing die associated wave- 
let bases of X^. Those wavelets provide a more efficient, 
non-redundant way of representing the difference images 
in Fig. 10. Since image reduction is achieved by repeated 
projection, the diflercncc between two successive signal 




J a Cubic spfine multiresolufion image approximation, (a) Cubic spfine pyramid {PoS,PiS ,P4S ,f^s}. 
(b) Error pyramid: {/Is - fts , ^iS - P45 , R, 5 - Hs. As }. What is displayed are the samples of those 
splines at the location of the knots (cardinal representation). 



approximations P^,-, / and 1^; f belong to the subspace 
W^- that is the complement of^S^" with respect to.Sj^., ; 

=5;, ew;]\vithS^; nW;\ ={0}.Thisisw1icrc 
the famous wavelct\j;(A') enters the scene: it generates the 
basis functions of the residual spaces [51], [91] 

T^;.=span{Hf(,v/2'-^)},,;,. 

There arc many applications (e.g., coding) where it is 
more concise to express die residues P^j., / " A' f^^^i' 
using wavelets rather rhan the basis functions ofV^".i as 
has been is done in Pig. 10. An example of wavelet trans- 
form is showai in Fig. 11; rliis decomposition works well 
for image coding because it produces many very small co- 
efficients in slowly var)ing image regions. 

In wavelet theory, splines constitute a case apart be- 
cause dicy give rise to the only wavelets that have a 
closed- form formula (pieccwise polynomial). All other 
wavelet bases arc defined indirccdy by an iiifmite re- 
cursion (or by an infinite product in the Fourier domain) 
[23], [48], [81], [109]. It is, therefore, no coincidence 
that most of the earlier wavelet constructiom were based 
on splines; tor instance, the Haar wavelet transform 
(«=0) [34], the Fi-anklin system (ff-1), Stromberg's 
one-sided orthogonal splines [82], and die celebrated 

Battlc-Lcmarid wavelets 
[11], [47]. Since then, the 
family has grown and there 
arc now several other sub- 
classes of spline wavelets 
available; diey differ in the 
type of projection used and 
in their orthogonality 
propcrdes. 

Corresponding to an or- 
tht>gonaI projection (and to 
thel^-pyi^amid above) is the 
class of scmi-OYtho£onal 
wavelets, which arc or- 
thogonal with rcspea to di- 
lation [98]. Tlicsc wavelets 
span the same space as tlic 
Batdc-Lcmarid splines, but 
ai*e not constrained to be or- 
Um^mml. This gives flexibil- 
ity and makes it ptissiblc to 
design wa\dets witii many 
interesting properties [5] 
and almost any desirable 
shape [1]. Of particular in- 
terest arc die B-spline wave- 
lets [20], [94], wliich arc 
compactly supported and 
optimally localized in time 
and frequency; asymptoti- 
cally, tiic}' adiievc the lower 
limit specified by Hciscn- 
bcrg's uncertaint)' principle 
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[94], The only downside of sani-oithogonal wavelets is diat 
some of the corresponding wavelet filters are UR. This is not 
a serious problem in practice, thanks to the availability of 
fast-recursive algoridims (cf. Box 2)—Kdealing efFiciendy with 
IIR filters is die main thinst of B-spline signal processing. 

Researchers have also designed spline wavelets such 
diat the con-esponding wavelet filters are FIR [21], [108], 
These biortho^onal wavelets are constructed using two 
multircsolutions instead of one, with the spline spaces on 
die syndiesis side. The major difference with the semi-or- 
thogonal case is diat die underlying projection operators 
are oblique rather than orthogonal [4]. Biortiiogonai 
spline \\^avelcts have many desirable properties that have 
made them very popular for applications: they are short, 
symmetrical, easy to implement (FIR filterbank), and very 
regular. Within the biordiogonal class, there is still one 
possibility'' which is to ortliogonalize die wavelets with re- 
spect to sliiftSj wloich leads to the more recent class of 
skfi-artho^oml wavelets. Such a construction was first il- 
lustrated with a femily of hybrid spline wavelets where the 
analysis and synthesis basis functions are spUnes of differ- 
ent degree andwj [104], 

Further Optimality Properties 
Variational Properties 

Splines have some very interesting extremal properties. 
One important result is the fint integral relation [3], 
which states that for any fiinaionf[x) whose mrh deriva- 
tive is square intcgrable, wc have 

-« _« -„ (35) 

where s{x) is the spline interpolaiit of degree «=2w-l 
such that =flk). In particular, if we apply this decom- 
position to the problem of the interpolation of a given 
data scqucnce/(^), we may conclude that, among all pos- 
sible intcrpolants/(:v;), the spline intcrpolant is the 
only one that minimizes the norm of the mt\\ deri\'ativc, 
which is a rather remarkable result [72], The reason is 
simply that the second term in (35) is non-zero if 
f{x):its{x) at the non-integer points. In this sense, die 
spline is the interpolating function that oscillates the least. 
For w=2, the cnci'gy function in (35) is a good approxi- 
mation to the integral of tiie curs^ature for the curve 
y=f(,x). Thus, cubic splLies intcrpolant exhibit a mini- 
mtm cmvmm property^ which justifies the analogy with 
the draftnian's spline, or French curve. The laaer device is 
a thin elastic beam that is constrained to pass through a 
given set of points, 

Smoothing Splines 

Interpolation is not the only approach for fitting a contin- 
uous model to a signal. For noisy data, an exact fit may 
not even be desirable. Such situations can be dealt with by 
relaxing the interpolation constraint and by maldng best 



use of our a priori laiowledge about the problem. The 
natural extension of the previous interpolation problem is 
to find the ftincrioni(A;) that minimizes 

X(/{*)-^W)^4-?/j(r'"'W)'^. 

klz L (36) 

This is a well-posed, regularized least-squares problem 
where the first term quantifies tiie error between the 
model s{x) and the measured data points/(^); tiie second 
term imposes a smootliness constraint on the solution. 
The choice of a particular value of the regularization fee- 
tor X reflects our a priori information; it can be based ei- 
ther on tiie Icnowledge of die variance of the noise or the 
degree of smootiiness of the signal as measured by (35). 
Here again, it can be shown that the optimal solution 
among all possible flinaions is a spline of degree n =2;»-l 
[65], [71]. Part ofthe argument follows from the first in- 
tegral equation: any non-spline fit can be improved by us- 
ing its spline interpolant that flirdier reduces the second 
term in the criterion while keeping the same y2.\uc&s{k) at 
the grid points. The solution to the above problem is 
called a mtoothing spline y because it is equivalent to a spe- 
cial form of smoothing of the data. Similar to the exact in- 
terpolation that corresponds to the case X— >0, the 
B-spline coefficients of die smoothing spline can be com- 
puted efficientiy by recursive filtering [96]. 

Introducing a regularization term, as in (36), is a stan- 
dard practice for dealing with many other types of 
ill-posed problems [61], including sparse and 
non-equally spaced data. The regularization parameter X 
is tj-pically used to control the smoothness of tiie solu- 
tion. For ?w=l, the regularization will tend to privilege 
small values of the derivative; a good physical analogy is 
that of a membrane that takes a constant value at cquilib- 




A / Separable cubic spline wavelet transform corresponding to 
tfie multiresolution decomposition in Fig. JO. The wavelet 
transform is implemented iterativeiy using a separable algo- 
rithm. First the rows ofthe image are split into two haifs using 
a two-channel filterbank. Second, the same procedure is ap- 
plied to the columns. This process is then iterated on the inter- 
mediate lower resolution images, whidi are precisely the ones 
displayed in Fig. JOa, 
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riutn. For w = 2, there is no pcnaJty for linear gradients. 
The generalization of this problem to higher dimensions 
leads to another area of study called "thin-plares splines" 
[113]. Generalized splines and mdiaJ basis functions can 
also be defined in a similar way by introducing more com- 
plex regularization terms [63], 

Smoothing splines are closely related to wavelet 
denoising techniques, which may be expressed in a regu- 
larization framework as well [19], The main difference is 
that the smoothing spline is a linear estimator, while 
Donoho's wavelet shrinkage is non-linear [30], The idea 
is simple and was pioneered by Weaver et al. using or- 
thogonal spline wavelets [114]; talce the wavelet trans- 
form of a signal and set to zero the coefficients below 
some critical threshold while slightly attenuating the 
odier ones (soft-threshold); tlien reconstruct the signal 
by inverse wavelet transform. The wavelet technique has 
the advantage of preser\ing edges; it is well suited for sig- 
nal or images that arc piecewise smooth, and is optimal in 
a wcU-defmed statistical sense [30], [48]. 

Best'ApproxImatlon Properties Among Wavelet 

In "Controlling the Approximation Error," we saw diat 
splines have mL-n-\-l order of approximation, which 
means that the error decays like theLth power of the sam- 
pling step. There are also non-spline fionctions ^>{x) that 
have the same property; in particular, thcLth order scal- 
ing funaions encountered in the multircsolution dieory 
of the wavelet transform. Note that, in the wavelet world, 
the order is usually specified by the number of vanishing 
moments of the anal)^is wavelet An equivalent 
statement of the order property is that the translates of the 
function <p must reproduce the polynomials of degi'ee n 
[26], [79] . In general, the order propert)^ implies that we 
have the following asymptotic fomi of the approximation 
error (cf. [89]) 

\\s - P,,r 4 = C ■ T ^ |, as r ^ 0, (37) 

where Pip.r^ is the projection of s onto the space 
Vt =sp^n{(^(x / T - k)} tez and where the constant C^,/, 
can be determined explicitly [14], [89], This is essentially 
the same equation as (24) with an equality instead of an 
upper bound; the. asymptotic leading constant C^^i is, 
therefore, necessarily smaller than Q in (24). 

Among ail known wavelet families, splines appear 
to have the best approximation property in the sense 
that the magnitude of the constant Ci^^i is minimum 
[83], [89]. This means that, in the asymptotic regime 
where the error is small, we can apply a coarser sam- 
pling step if we use splines as opposed to other basis 
functions (or wavelets) with the same order L, The po- 
tential reduaion in sampling density can be quite sig- 
nificant. For instance, Sweldens observed that splines 
at half the resolution could provide a better approxima- 
tion than Daubechies' wavelets [83]. Recently, the ex- 
act subsampling factor such that the asymptotic errors 



in both cases are identical has been determined analog- 
ically [14]; it converges to 7t as the order L gets suffi- 
ciendy large! 

Maximum Regularity and Shortest Support 

It is well known from wavelet theory that the B-splines 
are the shortest scaling flincdons of order Z [23], [81], 
They are also the most regular ones if one takes the size of 
die refinement filter into account [90] : their Sobolev reg- 
ularity (r derivatives in L2) is Vrrux =» + - [81] and their 
Holder exponent is a=» [66] . This latter property means 
that the B-spline of degree n is "almost" n times continu- 
ously-difFercntiabie; strictly speaking, the nth derivative 
of spline of degree n has some isolated points of disconti- 
nuities (knots), but is bounded nevertheless. 

If one extends the mathematical analysis to functions 
that do not necessarily satisfy the two-scale relation 
(multircsolution property), then the B-splines can still be 
shown to be the shortest ftinctions of order L. However, 
tiiere arc also other solutions, albeit less regular [12]. 
Thus, in the most general sense, the B-splincs are the 
shortest and smoothest functions of order I. Since the 
performance of an approximation algorithm is strongly 
determined by the order of approximation and to some 
extent by the regularity of the basis functions, this has im- 
portant practical consequences, especially for image in- 
terpolation (cf Box 3). In this type of processing, where 
computational cost is essentially determined by die size of 
the basis function, it makes pjerfect sense to use the short- 
est functions witli the required order properties; i.e., the 
B-splines. 



Fractional Splines 

Interestingly, B-splines can be generalized to fractional 
orders (cf. the illustration on the cover of this issue) 
[102], The fractional splines are piecewise power fimc- 
tions with building blocks of the form (x-Xi)°^ with 
a>— ^ real. The corresponding B-splines provide a 
smooth transition between the polynomial ones. They 
retain all the properties of the conventional 
B-splines — one merely replaces » by a in all formiilas, 
except the compact support [the finite sum in (10) be- 
comes an infinite one], One justification for looking at 
the fractional B-splines is that they oflfer the same con- 
ceptual ease for dealing with fractional derivatives as the 
conventional splines do for derivatives. One potential 
application is the analysis of fractional Brownian morion 
processes. 

Applications 

Our intent here is not to be exhaustive, but instead to give 
a brief overview of the type of signal and image process- 
ing applications that can benefit from the use of polyno- 
mial splines. 
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Zooming and Vlsoaltzatton 

Image zooming and interpolation are perhaps the most 
obvious applications of splines. These manipulations are 
especially useful for medical imaging [57], [60], but also 
for multimedia and digital photography, which arc rap- 
idly expandingapplications areas. The use of cubic splines 
in image processing was pioneered by Hou and Andrew 
[36], The proposed approach was not yet very practical 
because the B-spline coefficients were determined by ma- 
trix inversion. The method was made much more effi- 
cient with the introduction of recursive filtering 
algorithms [93]. Note that zooming by powers of two 
can also be implemented using the EXPAND function of 
a pyramid [96]. 

Geometric Image Transformations 

When there is no size reduction, geometric transforma- 
tions are often implemented by standard spline interpola- 
tion (cf Box 3). One of the drawbacks is that the 
complexity of tlie method, which is two-dimensional, 
grows rapidly with the order Z =n +1 of the model [typi- 
cally, 0{D) per pixel] . Fortunately, for the simplest trans- 
formations (scaling and rotations), there arc ways to 
make the problem separable through a clever 
factorization of die transformation matrix [58]. This 
technique was used in [105] to design a high-quality 
spline -based procedure allowing the rotation of images 
using ID convolutions only; it was extended in [86] to al- 
low for affine ti'ansformations in 2D and 3D as well. For 
image reductions, it is preferable to use a least-squares ap- 
proximation to reduce aliasing artifects. Such an algo- 
rithm exists for re-sizing images with arbitrar>' scaling 
factors [100] — not just the usual powers of two. Re- 
centiy, it has been simplified and accelerated using 
oblique projections [46]. The idea is to iise the box func- 
tion as the simplest possible prcfiltcr , and to apply the ap- 
propriate digital filtering compensation aftenyard so that 
die resulting approximation is a projcaion. The results 
are almost indistinguishable from the least-squares solu- 
tionj and the algorithm generalizes for any degree n. 

Filter Design and fast Continuous 
Wavelet Transform 

Thanlw to the wscale relation, a signal can be convolved 
very efficiendy with a disaete B-spline of size w using a 
cascade of moving average filters (recursive update). This 
jdelds an algorithm that has a complexity independent of 
the size of the basis functions. Thus, we have a very effi- 
cient way of implementing a scalable fil ter wiiose unpulsc 
response is the sum of a few B-spline basis functions. This 
is an idea that has been exploited for filter design [59], 
and for implementing the continuous wavelet transform 
with integer scales [101]. This type of algorithm achieves 
the lowest 0(1) complexity per corhputed coefficient. In 
contrast with other wavelet transform algorithms [67], 
the B-spline approach is non-iterative across scale and, 



therefore, well suited to a parallel implementation. 
Splines are also used to compute wavelet transforms with 
arbitrar)^ non-integer scales [110]. This is more compli- 
cated because it necessitates approximating enlarged 
wavelets using cither orthogonal [111] or oblique projec- 
tions [112]. This latter option appears to be more advan- 
tageous because it simplifies the determination of the 
filter coefficients without any measurable degi-adarion. 

Image Compression 

Image compression is anodier area where splines can be 
helpful. Most of today state-of-the-art methods use wave- 
lets — the most prominent ones are Shapiro's embedded 
zero-tree wavelet coder [77], and Said and Pearlman's 
SPIHT [69], While there are many possible choices of 
wavelet filters, many researchers tend to favor the 
biorthogonal splines for the reasons mentioned before 
(symmetry, short support, and excellent approximation 
properties) [9], [81]. We should also mention some 
non-wavelet-bascd systems; for example, the method of 
Toraichi et al., whidi uses quadratic spline interpolation 
[87], and Moulin's decomposition in terms of liierarchical 
spline basis functions [53]. Pyramid coders, which extend 
Burt and Adelson's initial idea, should not be dismissed ei- 
tiier [39], [56], [88], [107]. These can offer advantages, 
especially in higher dimensions where the overhead with 
respea to wavelets becomes negligible. Finally, splines 
provide a good solution for sub-pixel motion compensa- 
tion. Moulin et ai. have proposed a nicely integrated sys- 
tem where the motion vectore are represented using 
hierarchical basis functions (linear splines) [54]. 

MuHI'Scale Processing and Image Registration 

Spline pyramids provide a very convenient tool for per- 
forming mulriscale image proccssingj especially when tiie 
underlying problem is formulated in the continuous do- 
main. Tliis is a powerful idea for the implementation of it- 
erative algorithms using a coarse-to-fine iteration 
strategy [68], The benefits are twofold: first, there is an 
obvious acceleration because the cost of all 
low-resolution iterations is essentially negligible. Second, 
a multi-scale approach tends to be quite robustj wliich 
means that the algorithm is much less likely to get trapped 
in a local optimum. A good illustration of these ideas is 
provided by the image registration algorithm described 
in [85] . Tills method makes use of the same high-quality 
spline model for all aspects of the computation; image 
pyramid, geometric transform, and computation of the 
gradient of the criterion that is optimized. The benefits of 
this consistent design can be found in the results, w^hich 
are the best reported so far (error less than 1/lOOth of a 
pixel in a series of controlled experiments). The approach 
is reasonably fast because it makes the best use of its itera- 
tions: good starting conditions with an optimizer 
(Marquardt-Levenberg) that is extremely efficient near 
the optimum. 
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Contour Detection 

The spline formalism lends itself very naturally to the 
computation of gradients required for contour dctcaion. 
One can, for instance, reinterpret some of the classical 
edge detectors from this perspective [96] . To improve the 
gradient estimation in the presence of noise, Poggio et al. 
proposed using a smoothing spline technique [61], [62], 
They showed the approach to be more or less equivalent 
to smoothing the image widi a Gaussian filter in a prepro- 
cessing step (Cann/s edge detector [18]). This analogy 
holds even furdier [96]: there is an exact equivalence be- 
tween a smoothing spline edge detector and Deriche's re- 
cursive formulation of Canny's edge detector [28], 
Finally, Mallat and Zhong used wavelets that are deriva- 
tives of B -splines for obtaining their multi-scale edge rep- 
resentation of images [49]. 

Snakes and Contour Modeling 

In computer graphics^ curves are often generated using 
B-splines [10]. This parametric representation is also well 
suited to the analysis of shapes and contours [32], In par- 
ticular, it is well adapted to extracting shape invariants 
[22] , [37] . The simplest contour splines are piecewise lin- 
ear; they can be used to encode boundaries optimally in 
tlie rate-distorrion sense [44], [75], 

Mcnet et al. proposed using B-splincs snakes for ex- 
tracting contours in images [52]. A snake is an energy 
minimization spline segment vAth external and internal 
forces [43]. It simulates an elastic material that can dy- 
namically conform to local image features. The internal 
forces act as a regularizarion device by constraining the ri- 
gidity of the cm-ve. Alternatively, the smoothness of the 
curve can also be controlled direcdy and more simply by 
adapting the scale of the basis functions [16], 

Analog-tO'Dlgttal Conversion 

Spline and wavelet sampling present interesting alterna- 
tives to the conventional approach dictated by Shannon^s 
sampling theorem. These techniques can be adapted for 
dealing with non-ideal acquisition devices [92], and 
multi-channel measurements [106]. With tiiis more gen- 
eral view of samplmg, it is tempting to modify the acqui- 
sition scheme so as to measure the coefficients of some 
signal expansion (i.e., to perform some prescribed inner 
products) ratiier than to measure the samples of the signal 
itself Healy and Weaver have pioneered this idea for 
magnetic resonance imaging [35], [114]. They proposed 
a wavelet-encoding scheme using separable basis func- 
tions (Batde-Lemari^ splines along the A:-dimension, and 
conventional Fourier exponentials along die y-direction) . 
Splines are also useful for the converse task of digi- 
tai-to -analog conversion. Kamada et al. designed a qua- 
dratic spline signal generator [41], [42]; one of their 
circuits was used commercially for high-fidelity sound re- 
production. 



Conclusion 

We hope to have convinced the reader that splines consti- 
tute a useful tool for signal processing. Their main advan- 
tages can be summarized as follows; 
A One can always obtain a continuous representation of a 
discrete signal by fitting it with a spline in one or more di- 
mensions. The fit ma)^ be exact (interpolation) or approxi- 
mate {least-squares or smoothing splines). Spline fits are 
usually preferable to other forms of representations (e.g., 
Lagrange polynomial interpolation) because they have less 
of a tendency to oscillate (minimum curvature property) . 
A Polynomial splines can be expressed as linear combina- 
tions of B -spline basis functions. For equally spaced 
knots, the spline parameters (B -spline coefficients) may 
be determined by simple digital filtering. There is no need 
for matrix manipulations/ 

A The primary reason for worldng with the B -spline rep- 
resentation is that the B-splines are compacdy supported. 
They are the shortest functions with an order of approxi- 
mation!/ =« -I- 1. This short support property is a key con- 
sideration for computational efficiency. Their simple 
analytical form also greatly facilitates manipulations. 
A Splines are smooth and well-behaved functions 
(piecewise polynomials). Splines of degree n are («-l) 
continuously diflferentiable. As a result, splines have ex- 
cellent approximation properties. Precise convergence 
rates and error estimates are available. 
A Splines have mulriresolution properties that make 
them very suitable for constructing wavelet bases and for 
performing multi-scale processing, 
A B-splines and their wavelet coimterparts have excellent 
localization properties; they are good terhplates for 
time-frequency signal analysis, 

A The family of polynomial splines provides design flexi- 
bility. By increasing the degree », we progressively switch 
from the simplest . piecewise constant («=0) and 
piecewise linear (n-l) representations to the other ex- 
treme, which corresponds to a band-limited signal model 

(»^+oo). 

A The conventional sampling procedure can be easily 
modified to obtain a spline representation of an analog 
signal. This essentially amounts to replacing Shannon's 
ideal low-pass filter with another optimal prefiltcr speci- 
fied by the representation. In principlc> there is no com- 
pelling reason otiier. than liistory for preferring the 
band-limited model — and its corresponding sine 
interpolator — over other ones. 

Finally, similar spline techniques are also available for 
non-uniformly spaced data. The price to pay, however, is 
that one loses the convenient shift-invariant structure (fil- 
ters) that was emphasized in diis article. The reader who 
wishes to learn more about non- uniform splines is re- 
ferred to [24] and [74] . 
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